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Abstract. Finding a match between partially available deformable shapes 
is a challenging problem with numerous applications. The problem is 
O^l usually approached by computing local descriptors on a pair of shapes 

(— i and then establishing a point-wise correspondence between the two. In 

this paper, we introduce an alternative correspondence-less approach to 
matching fragments to an entire shape undergoing a non-rigid deforma- 
ns^ tion. We use diffusion geometric descriptors and optimize over the inte- 

gration domains on which the integral descriptors of the two parts match. 
The problem is regularized using the Mumford-Shah functional. We show 
an efficient discretization based on the Ambrosio-Tortorelli approxima- 
tion generalized to triangular meshes. Experiments demonstrating the 
success of the proposed method are presented. 

Keywords: deformable shapes, partial matching, partial correspondence, 
1 1 partial similarity, diffusion geometry, Laplace-Beltrami operator, shape 

descriptors, heat kernel signature, Mumford-Shah regularization 
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00 1 Introduction 

wn In many real-world settings of the shape recognition problem, the data are de- 

y—( graded by acquisition imperfections and noise, resulting in the need to find partial 

similarity of objects. Such cases are common, for example, in face recognition, 
where the facial surface may be partially occluded by hair. In other applications, 
' . I such as shape retrieval, correct semantic similarity of two objects is based on 

^ partial similarity - for example, a centaur is partially similar to a human because 

k>( they share the human- like upper body |15j . 

j_j In rigid shape analysis, modifications of the popular iterative closest point 

(ICP) algorithm are able to deal with partial shape alignment by rejecting points 
with bad correspondences (e.g., by thresholding the product of local normal vec- 
tors). However, it is impossible to guarantee how large and regular the resulting 
corresponding parts will be. 

Bronstein et al. [4 formulated non-rigid partial similarity as a multi-criterion 
optimization problem, in which one tries to find the corresponding parts in two 



2 



Pokrass et al. 



shapes by simultaneously maximizing significance and similarity criteria (in [4], 
metric distortion |14|19l5j was used as a criterion of similarity, and part area 
as significance) . The problem requires the knowledge of correspondence between 
the shapes, and in the absence of a given correspondence, can be solved by 
alternating between weighted correspondence finding and maximization of part 
area. In [4;, a different significance criterion based on statistical occurrence of 
local shape descriptors was used. 

One of the drawbacks of the above method is its tendency in some cases to 
find a large number of disconnected components, which have the same area as 
a larger single component. The same authors addressed this problem using a 
Mumford-Shah |22ll0j -like regularization for rigid [3J and non-rigid [5] shapes. 

Recent works on local shape descriptors (see, e.g., descriptors 16 33 [24|11)18|2 1 32 29 27|9j ) 
have led to the adoption of bags of features [28] approach popular in image anal- 
ysis for the description of 3D shapes |21I23I30| . Bags of features allows to some 
extent finding partial similarity, if the overlap between the parts is sufficiently 
large. 

In this paper, we present an approach for correspondence-less partial match- 
ing of non-rigid 3D shapes. Our work is inspired by the recent work on partial 
matching of images |12j . The main idea of this approach, adopted here, is to find 
similar parts by comparing part-wise distributions of local descriptors. This re- 
moves the need of correspondence knowledge and greatly simplifies the problem. 

The rest of the paper is organized is as follows. In Section 2, we review the 
mathematical background of diffusion geometry, which is used for the construc- 
tion of local descriptors. Section 3 deals with the partial matching problem and 
Section 4 addresses its discretization. Section 5 presents experimental results. 
Finally, Section 6 concludes the paper. 

2 Background 

Diffusion geometry. Diffusion geometry is an umbrella term referring to ge- 
ometric analysis of diffusion or random walk processes. We models a shape as 
a compact two-dimensional Riemannian manifold X. In it simplest setting, a 
diffusion process on X is described by the partial differential equation 



called the heat equation, where A denotes the positive-semidefinite Laplace- 
Beltrami operator associated with the Riemannian metric of X. The heat equa- 
tion describes the propagation of heat on the surface and its solution f(t,x) is 
the heat distribution at a point x in time t. The initial condition of the equa- 
tion is some initial heat distribution f(0,x); if X has a boundary, appropriate 
boundary conditions must be added. 

The solution of corresponding to a point initial condition /(0, x) = S(x, y) : 
is called the heat kernel and represents the amount of heat transferred from x 
to y in time t due to the diffusion process. The value of the heat kernel ht(x,y) 




(1) 
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can also be interpreted as the transition probability density of a random walk 
of length t from the point x to the point y. 

Using spectral decomposition, the heat kernel can be represented as 



Here, fa and A^ denote, respectively, the eigenfunctions and eigenvalues of the 
Laplace-Beltrami operator satisfying A fa — \i<fii (without loss of generality, 
we assume Xi to be sorted in increasing order starting with Ao = 0). Since 
the Laplace-Beltrami operator is an intrinsic geometric quantity, i.e., it can be 
expressed solely in terms of the metric of X, its eigenfunctions and eigenvalues as 
well as the heat kernel are invariant under isometric transformations (bending) 
of the shape. 

Heat kernel signatures. By setting y — x, the heat kernel h t (x,x) ex- 
presses the probability density of remaining at a point x after time t. The value 
ht(x, x), sometimes referred to as the auto-diffusivity function, is related to the 
Gaussian curvature K(x) through 



This relation coincides with the well-known fact that heat tends to diffuse slower 
at points with positive curvature, and faster at points with negative curvature. 
Under mild technical conditions, the set {h t (x, a;)}t>o is fully informative in the 
sense that it allows to reconstruct the Riemannian metric of the manifold [29] . 

Sun et al. |29] proposed constructing point-wise descriptors referred to as 
heat kernel signatures (HKS) by taking the values of the discrete auto-diffusivity 
function at point x at multiple times, p(x) = c{x){h tl (x, x), ...,ht d (x,x)), where 
ti, td are some fixed time values and c(x) is chosen so that ||p(a;)||2 = 1 • Such 
a descriptor is a vector of dimensionality d at each point. Since the heat kernel 
is an intrinsic quantity, the HKS is invariant to isometric transformations of the 
shape. 

A scale- invariant version of the HKS descriptor (SI-HKS) was proposed in [9]. 
First, the heat kernel is sampled logarithmically in time. Next, the logarithm and 
a derivative with respect to time of the heat kernel values are taken to undo the 
multiplicative constant. Finally, taking the magnitude of the Fourier transform 
allows to undo the scaling of the time variable. 

Bags of features. Ovsjanikov et al. [23] and Toldo et al. [30] proposed con- 
structing global shape descriptors from local descriptors using the bag of features 
paradigm [28]. In this approach, a fixed "geometric vocabulary" is computed by 
means of an off-line clustering of the descriptor space. Next, each point de- 
scriptor is represented in the vocabulary using vector quantization. The bag of 
features global shape descriptor is then computed as the histogram of quantized 
descriptors over the entire shape. 




(2) 



i>0 
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3 Partial matching 

In what follows, we assume to be given two shapes X and Y with correspond- 
ing point-wise descriptor fields p and q defined on them (here we adopt HKS 
descriptors, though their quantized variants or any other intrinsic point-wise de- 
scriptors can be used as well). Assuming that Y is a part of an unknown shape 
that is intrinsically similar to X, we aim at finding a part X' C X having the 
same area A of Y such that the integral shape descriptors computed on X' and 
Y coincide as closely as possible. In order to prevent the parts from being frag- 
mented and irregular, we penalize for their boundary length. The entire problem 
can be expressed as minimization of the following energy functional 



E(X') = 



pda — q 



X' 



+ X T L(dX') 



(4) 



under the constraint A(X') = A, where A denotes area and q = / qda. The 

first term of the functional constitutes the data term while the second one is the 
regularity term whose influence is controlled by the parameter A r . 

Discretization of the above minimization problem with a crisp set X' re- 
sults in combinatorial complexity. To circumvent this difficulty, in |2|3j it was 
proposed to relax the problem by replacing the crisp part X' by a fuzzy mem- 
bership function aonl, replacing the functional E by a generalization of the 
Mumford-Shah functional [22] to surfaces. Here, we adopt this relaxation as well 
as the approximation of the Mumford-Shah functional proposed by Ambrosio 
and Tortorelli [1] . This yields the problem of the form 



min D(u) + A T R(u; p) 

u,p,a 

s.t. / uda = A, 
Jx 



(5) 



with the data term 



D(u) = 



puda — q 



x 



and the Ambrosio- Tortorelli regularity term 
A 



(6) 



R(u;p) 



I p 2 \\Vu\\ 2 da + \ h e [ ||Vp|| 2 da+^ / (1 - pfda, (7) 
Jx Jx 4e J x 



where p is the so-called phase field indicating the discontinuities of u, and e > 
is a parameter. 

The first term of R above imposes piece-wise smoothness of the fuzzy part u 
governed by the parameter A s . By setting a sufficiently large A s , the parts become 
approximately piece-wise constant as desired in the original crisp formulation 
Q. The second term of R is analogous to the boundary length term in Q and 
converges to the latter as e —} 0. 
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We minimize §5§ using alternating minimization comprising the following two 
iteratively repeated steps: 
Step 1: fix p and solve for u 



puda — q 



x 



+ A r 



p 2 ||Vu|| 2 da s.t. / uda = A 



x 



x 



Step 2: fix the part u and solve for p 
. A 



mm — 
p 2 



x 



A,. 



p z \\Vu\\ z da + A b e / \\V p\\ 2 da + / (f - pfda 



x 



4c 



(8) 



(9) 



x 



4 Discretization and numerical aspects 

We represent the surface X as triangular mesh with n faces constructed upon the 
samples {xi, . . . , x m } and denote by a = (ai, . . . , a m ) T the corresponding area 
elements at each vertex. A = diag{a} denote the diagonal to x to matrix created 
out of a. The membership function u is sampled at each vertex and represented 
as the vector u = {u\, . . . ,u m ) T . Similarly, the phase field is represented as the 
vector p = (/?!,.. . ,p m ) T . 

Descriptors. The computation of the discrete heat kernel /i t (xi,X2) re- 
quires computing discrete eigenvalues and eigenfunctions of the discrete Laplace- 
Beltrami operator. The latter can be computed directly using the finite elements 
method (FEM) [26], of by discretization of the Laplace operator on the mesh 
followed by its eigendecomposition. Here, we adopt the second approach accord- 
ing to which the discrete Laplace-Beltrami operator is expressed in the following 
generic form, 

(Af)i = ^I>i(/i-/i), (io) 

where /j = /(x^) is a scalar function defined on the mesh, are weights, 



and at are normalization coefficients. In matrix notation, ( 10 ) can be written as 
At = A _1 Wf , where f is an to x I vector and W = diag w f| — w ij- 

The discrete eigenfunctions and eigenvalues are found by solving the gener- 
alized eigendecomposition [17 W<& = A3? A, where A = diag{A;} is a diagonal 
matrix of eigenvalues and $ — (<pi(xi)) is the matrix of the corresponding eigen- 
vectors. 

Different choices of W have been studied, depending on which continuous 
properties of the Laplace-Beltrami operator one wishes to preserve [13131] . For 
triangular meshes, a popular choice adopted in this paper is the cotangent weight 
scheme [25120] . in which 

w = f (cotAi + cot7y)/2: x J eAA(x i ); 
lJ \ : else, 

where and jij are the two angles opposite to the edge between vertices x^ 
and Xj in the two triangles sharing the edge. 
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Data term. Denoting by P d x m the matrix of point-wise descriptors on 
X (stored in columns), we have 



puda w Pdiag{A}u 



This yields the following discretization of the data term Q 



D(u) = ||PAu- qf 



m rn m rri m 

u A P PAu - 2q T PAu + q T q. 



(12) 



(13) 



Gradient norm. We start by deriving the discretization of a single term 
/9 2 ||Vw|| 2 da in some triangle t of the mesh. Let us denote by Xi,Xj and the 
vertices of the triangle and let X( = (xj — Xi,Xfc — x^) be the 3x2 matrix 

whose columns are the vectors forming the triangle, and by a t — \ \J det(X^X t ) 
its area. Let also D t be the sparse 2 x m matrix with +1 at indices and 
(2,fc), and —1 at (l,i) and (2,i). D t is constructed in such a way to give the 
differences of values of u on the vertices of the triangle with respect to the 
values at the central vertex, D t u = (uj — ui.Uk — Ui) T . The gradient of the 
function u is constant on the triangle and can be expressed in these terms by 
St = (Xj'XtJ-^DtU = G t u. 

In order to introduce the weighting by p 2 , let S be an n x m sparse matrix 
with the elements su — \ for every vertex i belonging to the triangle t and 
zero otherwise. In this notation, Sp is a per-triangle field whose elements are 
the average values of p on each of the mesh triangles. We use the Kroenecker 
product of S with 1 = (1, 1) T to define the 2n x n matrix S ® 1 formed by 
replicating twice each of the rows of S. This yields 



J p 2 ||V U || 2 da«]T||(Sp) t G t u|| 



diag{(S ® l)p} 



at 

2 



u T G T S(p)Gu, 



(14) 



where G is the matrix containing ^JalGt stacked as rows, and S(p) = diag{(S<8) 
1)M 2 - 

Discretized alternating minimization. We plug in the results obtained 
so far into the two steps of the alternating minimization problem ([8|~([9]). For 
fixed p, the discretized minimization problem Q w.r.t. u can be written as 

minu T ^A T P T PA + A r yG T S(p)G^ u-2q T PAu s.t. a T u = A(15) 

Let us now fix u. In a triangle t, we denote g 2 = ||Gtu|| 2 and let R(u) = 
diag{ai<? 2 , . . . , a n g 2 }. Using this notation, we obtain the following discretization 
of the integrals in the regularization term Q 



f P 2 \\\7u\\ 2 da « ^(Sp)?^a ( = p T S T R(u)Sp. 
Jx t 



(16) 
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Similar to the derivation of ( 14 ) , 



J ||V i o|| 2 rfawp T SG T GSp, (17) 



and 



/ (1 - pfda w p T Ap - 2a T p + l T a. (18) 
Jx 

The discretized minimization problem ^ w.r.t. p becomes 

min p T (2eA s S T R(u)S + 4e 2 A b SG T GS + A b A) p - 2A b a T p. (19) 



Since the above is an unconstrained quadratic problem, it has the following 
closed-form solution 

p= ( 2 ^S T R(u)S + 4 e 2 SG T GS + A) 'a. (20) 



5 Results 

In order to test our approach, we performed several partial matching experiments 
on data from the SHREC 2010 benchmark 07] and the TOSCA dataset [5]QThe 
datasets contained high-resolution (10K-50K vertices) triangular meshes of hu- 
man and animal shapes with simulated transformations, with known groundtruth 
correspondence between the transformed shapes. In our experiments, all the 
shapes were downsampled to approximately 2500 vertices. Parts were cut by 
taking a geodesic circle of random radius around a random center point. 

For each part, the normalized HKS descriptor was calculated at each vertex 
belonging to the part. To avoid boundary effects (see Figure [TJ, descriptors close 
to the boundary were ignored when calculating q in Q. The distance from the 
boundary was selected in accordance to the time scales of the HKS. We used ten 
linearly spread samples in range [65,90] for the descriptors and the according 
distance taken from edge was set to 15. Two to three iterations of the alternating 
minimization procedure were used, exhibiting fast convergence (Figure [2]). After 
three iterations the member function u typically ceased changing significantly. 
The phase map p assumed the values close to 1 in places of low gradient of the 
membership function u, and less than 1 in high gradient areas (Figure |J). The 
importance of the regularization step is is evident observing the change in u 
in Figure [2j Figure [3] shows the influence of the parameter A r , controlling the 
impact of the regularization. For too small values of A r , two equally weighted 
matches are obtained due to symmetry (left). The phenomenon decreases with 
the increase of the influence of the regularization penalty. However, increasing A r 
more causes incorrect matching (second and third columns from the right) due 



1 Both datasets are available online at http://tosca.cs.technion.ac.il 
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to low data term influence. Increasing it even more starts smoothing the result 
(rightmost column) until eventually making the membership function uniform 
over the entire shape. The resulting membership function u was thresholded in 
such a way that the outcome area will be as close as possible to the query area. 
Figures [4}j6] show examples of matching results after thresholding. Notice that 
in FigureH] the matching result sometimes contain the symmetric counter part 
of the result due to invariance of the HKS descriptor to intrinsic symmetry (in 
this figure, the threshold was adjusted to the value of 0.35 when the membership 
functions weights are split between two symmetric parts as in Figure [3|. The 
method is robust to shape deformations and geometric and topological noise as 
depicted in Figure [6j Note that the figures show part-to-whole shape match- 
ing, but because of the low scale HKS descriptors the same procedure works 
for matching to other parts as well. Table [l] summarizes quantitative evaluation 












- 







Fig. 1. An RGB visualization of the first three component of HKS descriptors com- 
puted on the full shape (left) and on a part of the shape (center). The L2 difference 
between the two fields is depicted in the rightmost figure. Note that the difference is 
maximal on the boundary decaying fast away from it; the error decay speed depends 
on the scale choice in the HKS descriptor. 



that was performed on a subset of the SHREC database, for which groundtruth 
correspondence and its bilaterally symmetric counterpart were available. This 
subset included a male, a dog and a horse shape classes with different geomet- 
ric, topological and noise deformations (98 shapes in total). The query set was 
generated by selecting a part from a deformed shape (1000 queries in each defor- 
mation category) and matched to the null shape with parameters and thresholds 
as described above. 

Complexity. The code was implemented in Matalb with some parts written in 
C with MEX interface. The quadratic programming problem ^ in Step 1 was 
solved using QPCj^Jimplementation of a dual active set method. The experiments 
were run on 2.3GHz Intel Core2 Quad CPU, 2GB RAM in Win7 32bit environ- 



2 available online at 



http : / / sigpromu . org/quadprog 
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Fig. 3. The influence of the parameter A r , controlling the impact of regularization. The 
leftmost figure depicts a query part; figures on the right are the membership function 
u for different values of A r . 



merit. The running time (including re-meshing and descriptor calculation) per 
part was 40 — 50 sec. 



6 Conclusions 



We presented a framework for finding partial similarity between shapes which 
does not rely on explicit correspondence. The method is based on regularized 
matching of region- wise local descriptors, and can be efficiently implemented. Ex- 
perimental results show that our approach performs well in challenging match- 
ing scenarios, such as the presence of geometric and topological noise. In the 
future work, we will extend the method to the setting of two partially-similar 
full shapes, in which two similar parts have to be found in each shape, and then 
consider a multi-part matching (puzzle) scenario. 
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Table 1. Part matching performance on transformed shapes from the SHREC bench- 
mark. At each query a random part (location and size) was selected from a deformed 
shape and matched to the null shape. Overlap is reported compared to the groundtruth 
correspondence between the shapes (in parentheses taking into account the intrinsic 
bilateral symmetry). 



Transformation 


Queries 


Avg. overlap 


Isometry 


1000 


75% (85%) 


Isometry + Shotnoise 


1000 


75% (85%) 


Isometry + Noise 


1000 


71% (82%) 


Isometry + Microholes 


1000 


68% (82%) 


Isometry + Holes 


1000 


66% (76%) 


All 


5000 


71% (82%) 
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Fig. 4. Examples of matching of random parts of shapes (first row) to approximately 
isometric deformations of the shapes (second row). Color code indicates different parts. 




Fig. 5. Examples of matching of random parts of shapes (first row) to to approximately 
isometric deformations of the shapes (second row). 




Fig. 6. Results of matching a part of shape (first row) to shapes distorted by different 
transformations (second row). Shown left-to-right are: shot noise, noise, micro holes 
and holes. 



